/*
 * Main.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include "Applications/Equilibrator.h"
#include "Multiphase.h"
#include "ReactionSystem.h"
#include "EquilibriumSolver/EquilibriumSolver.h"
#include "KineticSolver.h"

#include <iomanip>
#include "Utils.h"
#include "Test/Activity/HkfModelTest.h"
#include "Test/TRISpH/TRISpH.h"
#include "Utilities/Range.h"
#include "Numerics/LinearInterpolator.h"
#include "Test/Numerics/LagrangeInterpolatorTest.h"
#include "Test/Numerics/LinearInterpolatorTest.h"
#include "Assembler.h"
#include "System.h"
#include "Applications/Reactor.h"
#include "Demo/HaliteSolubility.h"
#include "Demo/CO2Solubility.h"

void TotalEquilibriumSolverTest();

int main()
{
	//DemoHaliteSolubility();
	
	DemoCO2Solubility();
	
	//TotalEquilibriumSolverTest();
	
//	System system = AssembleSystem({"H2O", "NaCl(aq)", "Calcite"});
//	
//	Equilibrator equilibrator;
//	
//	equilibrator.AddCondition("ChargeBalanced");
//	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
//	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
//	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
//	equilibrator.AddCondition("FixedTotalMass", "Ca[2+]", 2.0);
//	equilibrator.AddCondition("FixedTotalMass", "HCO3[-]", 2.0);
//	//equilibrator.AddCondition("pH", 10.060582);
//	
//	equilibrator.Equilibrate(system);
//	
//	cout << system << endl;
	
//	System system = AssembleSystem({"H2O", "NaCl(aq)", "Calcite"}, {"Calcite"});
//	
//	system.SetComposition("Calcite", 2.0);
//	
//	Equilibrator equilibrator;
//	
//	equilibrator.AddCondition("ChargeBalanced");
//	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
//	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
//	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
//	equilibrator.AddCondition("FixedTotalMass", "Ca[2+]", 1.0E-7);
//	equilibrator.AddCondition("FixedTotalMass", "HCO3[-]", 1.0E-7);
//	
//	equilibrator.Equilibrate(system);
//	
//	cout << system << endl;
//	
//	Reactor reactor(system);
//	
//	double t = 0.0, tEnd = 1000.0;
//	
//	cout << setprecision(6) << fixed;
//
//	while(t < tEnd)
//	{
//		cout << setw(15) << left << t;
//		cout << setw(15) << left << system.Get("Moles", "Calcite");
//		cout << setw(15) << left << system.Get("KineticRate", "Calcite");
//		cout << setw(15) << left << system.Get("KineticRate", "Ca[2+]");
//		cout << setw(15) << left << system.Get("pH") << endl;
//		
//		reactor.React(t, tEnd);
//	}
	
	
//	double tEnd = 1000.0, dt = 1.0E-3, t = 0.0;
//
//	  double T = system.GetTemperature();
//	  double P = system.GetPressure();
//	VectorXd n = system.GetComposition();
//	
//	KineticSolver kSolver(system.GetMultiphase(), system.GetEquilibriumReactions(), system.GetKineticReactions());
//	
//	kSolver.SetControlSolution(1.0E-6);
//	
//	cout << setprecision(6) << fixed;
//	
//	while(t < tEnd)
//	{
//		cout << setw(15) << left << t;
//		cout << setw(15) << left << system.Get("Moles", "Calcite");
//		cout << setw(15) << left << system.Get("KineticRate", "Calcite");
//		cout << setw(15) << left << system.Get("KineticRate", "Ca[2+]") << endl;
//		
//		kSolver.Iterate(t, tEnd, dt, T, P, n);
//		
//		system.SetComposition(n);
//	}
	
	//LinearInterpolatorTest();
	
	//LagrangeInterpolatorTest();
	
	//TRISpH(25.0, 1.0);
	
	//TRISpHTRange(Range(25.0, 150.0, 1.0), 1.0, 1.0, 0.01);
	
	//NoTRISpHTRange(Range(25.0, 150.0, 1.0), 1.0, 1.0E-7);
	
	//TestHkfModel_WaterNaCl(25.0, 1.0);
	
	//TestHkfModel_WaterCaCl2(25.0, 1.0);
}

void TotalEquilibriumSolverTest()
{
	System system = AssembleSystem({"H2O", "NaCl(aq)"});
	
	Equilibrator equilibrator;
	
	equilibrator.AddCondition("ChargeBalanced");
	equilibrator.AddCondition("FixedTotalMass", "H2O", 55.508);
	equilibrator.AddCondition("FixedTotalMass", "Na[+]", 1.0);
	equilibrator.AddCondition("FixedTotalMass", "Cl[-]", 1.0);
	
	equilibrator.Equilibrate(system);
	
	cout << system << endl;
	
	const Multiphase& multiphase = system.GetMultiphase();
	
	const CanonicalReactionSystem& reactions = system.GetEquilibriumReactions();
	
	cout << reactions.GetPrimarySpecies() << endl;
	
	TotalEquilibriumSolver totalSolver(multiphase, reactions);
	
	VectorXd n = VectorXd::Constant(multiphase.GetSpecies().size(), 1.0E-7);
	
	n[system["H2O"]] = 35.508;
	n[system["Na[+]"]] = 2.0;
	n[system["Cl[-]"]] = 2.0;
	
	VectorXd uj(4);
	
	uj << 1.0, 55.508, -1.0E-7, 1.0;
	
	totalSolver.Solve(uj, 25.0, 1.0, n);
	
	system.SetComposition(n);
	
	cout << system;
}
